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Abstract 

An algorithm is presented, that provides a fast and robust reconstruction of neutrino induced upward-going muons 
and a discrimination of these events from downward-going atmospheric muon background in data collected by the 
ANTARES neutrino telescope. The algorithm consists of a hit merging and hit selection procedure followed by fitting 
steps for a track hypothesis and a point-like light source. It is particularly well-suited for real time applications such as 
online monitoring and fast triggering of optical follow-up observations for multi-messenger studies. The performance 
of the algorithm is evaluated with Monte Carlo simulations and various distributions are compared with that obtained 
in ANTARES data. 



1. Introduction 

The main goal of neutrino telescope experiments such as AMANDA 111], IceCube iQl, NT-200 in lake Baikal Jst] 
and ANTARES is the observation of high energy neutrinos from non-terrestrial sources. These instruments detect 
Cherenkov light from the passage of relativistic charged particles produced in neutrino interactions in the detector or 
its surrounding material. By measuring the photon arrival times at known positions, the particle trajectory can be re- 
constructed. The flight direction of these particles is nearly colinear with the incident neutrino direction for neutrinos 
above 10 TeV. The best angular resolution can be reached for v^, charged current interactions where the measured 
particle is a muon that can travel several kilometres in water at TeV energies. Events induced by neutrinos from 
astrophysical sources must be distinguished from background which originates in the Earth's atmosphere. Whereas 
atmospheric neutrinos are considered a non-reducible background (at least without the use of their energy), atmo- 
spheric muons can in principle be excluded by simple geometrical considerations. To exclude the contribution from 
downward-going atmospheric muons, it is sufficient to identify upward-going events, which can only be produced 
by neutrinos. This distinction requires a reliable determination of the elevation angle, because downward-going at- 
mospheric muons outnumber upward-going atmospheric neutrinos by 5 to 6 orders of magnitude for typical neutrino 
telescope installation depths. 

In the following, a tracking algorithm is presented which meets this goal. Apart from providing an excellent up- 
down separation, it is very fast and therefore well-qualified for real time applications. This has been demonstrated in 
ANTARES, where the algorithm can cope with a trigger rate of 1 00 Hz, running on a single PC. It is used in various 
online monitoring tasks, for an analysis of atmospheric muons |@[| and as an alert sending program to trigger optical 
follow-ups of selected neutrino events Other potential applications include the study of magnetic monopoles 
and nuclearites and the measurement of atmospheric neutrinos. Thepresented algorithm is complementary to other 
reconstruction methods which have been developped in ANTARES ijx^. 

Section 2 describes the ANTARES detector. Section 3 explains the geometrical approximations which are used in 
all subsequent steps. The hit merging and hit selection methods are detailed in Sections 4 and 5 whereas the fitting 
procedure can be found in Section 6. The performance of the entire algorithm introduced in Sections 2-6 is discussed 
in Section 7. 

2. The ANTARES detector 

The ANTARES detector is located about 40 km off the coast of Toulon, France, at a depth of 2475 m in the 
Mediterranean Sea (42°48'N, 6°10'E). It consists of 12 flexible lines, each with a total height of 450 m, separated by 
distances of 60-70 m from each other They are anchored to the sea bed and kept near vertical by buoys at the top 
of the lines. Each line carries a total of 75 10-inch Hamamatsu photo multipliers (PMTs) housed in glass spheres, 
the optical modules (OM) |@], arranged in 25 storeys (three optical modules per storey) separated by 14.5 m, starting 
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100 m above the sea floor. Each PMT is oriented 45° downward with respect to the vertical. A titanium cylinder in 
each storey houses the electronics for readout and control. 

The positions of the line anchors on the sea bed and the construction details (cable length etc.) of each detector line 
are well known and stable in time. However the absolute positions and orientation of the individual OMs vary with 
time as the lines move freely in the sea current. Their movements are monitored by a system of acoustic transponders 
and receivers distributed over the lines and on the sea bed. After each cycle of acoustic signal exchange between 
these elements, a triangulation provides the shape of each line. This is complemented by data from a compass card 
within each electronics cylinder which measures the heading and the inclination of its storey. The processing of these 
combined data provides the position and orientation of each OM with a precision of better than 10 cm, but after an 
important delay. 

3. Geometrical approximations 

To qualify for real time applications the presented method uses only time independent geometrical informations 
and the following approximations are applied: 

Line shape Detector lines are considered to be perfectly vertical. Line distortions due to sea currents are ignored. 

Storey geometry The detailed geometry of the storey is ignored. Each storey is considered as a single OM, which is 
located directly on the detector line with an axis-symmetric field of view. The signals of the three PMTs within 
one storey are combined. 

These approximations are used in the hit merging, hit selection and fitting steps. Their effect on the angular resolution 
is discussed in Section [73] Under the above approximations, a detector line is a straight line in space. It is also 
assumed that a muon track is a straight line, i.e. multiple scattering etc. is ignored. Except for the special case that 
the muon track and a detector line are parallel (in the case of an exactly vertical track), the point of closest approach 
of the muon track to the detector line can be determined. Most of the Cherenkov light is expected to be seen in the 
vicinity of this point. This is used in the hit selection as well as in the fitting procedure of the algorithm described 
below. 



4. Hits and hit merging 

The PMT signals are digitized by a custom built ASIC chip ifioll . If the analog PMT signal crosses a preset 
threshold, typically 1/3 photoelectron (p.e.), its arrival time is measured together with its charge. The latter is obtained 
by integrating the analog PMT signal within a time window of 40 ns. Each such pair of time and charge is called a 
hit and the corresponding data for each hit are sent to shore. The data stream is processed by a computer farm in the 
shore station which searches for different physics signals according to predefined trigger conditions. If at least one 
trigger pattern is found, all hits within a time window of few //s around the trigger time are stored on disk. The DAQ 



system is described in detail in 111 111 . 

For the data which will be presented in Section|7]3]two trigger algorithms are used. Both rely on LI hits which are 
defined either as 2 hits in coincidence within 20 ns in two OMs on the same storey or as a single hit with an amplitude 
larger than a threshold of 3 or 10 p.e.. One of the used triggers requires at least five causally connected LI hits with 
respect to a given track direction. The second used trigger requires two pairs of LI hits in adjacent or next-to-adjacent 
storeys within 80 ns and 160 ns respectively. 

Before starting the actual hit selection, time and charge calibrations are applied to transform raw data to physical 
units such as time (in nanoseconds) and charge (in equivalent number of photoelectrons). Up to now the standard 
data handling in ANTARES has been described. The following steps are instead related to the presented algorithms. 
Different reconstruction methods might apply other hit treatment and hit selection procedures. 

The calibrated hits are time-ordered for each storey. All hits on the same storey within a predefined time window 
are merged, as applied to early ANTARES data in 0]. Merging is performed by adding the hit charges and keeping 
the time of the earliest hit. The merging time window has to be chosen short enough to ensure a clean hit selection, 
with the earliest of the merged hits having a high probability of originating from a particle track. On the other hand. 
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the merging time window must be long enough to select most of the signal hits. For a 10-inch phototube, operated 
at a threshold of 0.3 photoelectrons a typical noise rate of 60-100 kHz is observed in sea water lll2ll . In this case a 
merging window of 20 ns is appropriate. As a result, a list of merged hits per storey is obtained and is the basis for all 
subsequent steps. 



5. Hit selection 

The light pattern expected from a highly relativistic charged particle is illustrated in Figure [1] The goal of the hit 
selection is to choose hits due to Cherenkov photons and avoid random hits from the optical background or scattered 
late hits. Late hits due to scattering are rare compared to typical optical background of 60-100 kHz. The probability 
of a Cherenkov photon at A - 470 nm to suffer from light scattering on a path of 60 m (the inter-line distance and 
coincidentally the water absorption length) is only few per cent ifisll . 

For tracks which are not parallel to a detector line, most of the Cherenkov Ught is expected near the point of closest 
approach of the track to the detector line. Therefore, the first step of the hit selection consists of finding a "hot spot" 
of light on a detector line which originates from Cherenkov light induced by a passing particle. In the presence of 
background light, a single high charge hit does not provide a "hot spot" with sufficient purity; two high charge hits 
in adjacent or next-to-adjacent storeys are needed. The detector storeys are numbered consecutively along the z-axis, 
like floors in a building. For a given floor ; the adjacent or next-to-adjacent storeys are at floor ; + j with j = 1, 2. To 
ensure that hits on floors ; and / + j were caused by the same particle, the following condition on the absolute hit time 
difference At is imposed 

Af < jAz- + t„ (1) 
c 

with Az the absolute vertical distance between adjacent storeys, c the speed of light in vacuum, « the refractive index 
of the medium and f , an additional time delay which accounts for the timing uncertainty of the photon detection 
itself and the simplifications of the detector geometry. It is assumed that storeys are equally spaced. Such a causality 
condition has been used for hit selection prior to track reconstruction in Baikal ^ and AMANDA lfl4ll . For the 
ANTARES detector with Az - 14.5 m and fj = 10 ns, one finds Af < 80 ns for adjacent and At < 150 ns for next- 
to-adjacent storeys. Only merged hits which originate from different PMTs or with a charge above 2.5 photoelectrons 
are considered here, because they have a high probability of being caused by a particle track. All hits which contribute 
to the above defined "hot spots" are stored as selected hits for the next steps. If several "hot spots" occur in the same 
storey at different times, only the earliest of them is kept. 

For the following steps, only detector lines having at least one "hot spot" are considered. This condition removes 
all isolated hits from the sample. Using the "hot spot" hits as seeds, hits are added to the list of selected hits. From 
this point on, the hit selection is exclusively based on timing information. No charge cut is applied at this stage. This 
allows to use the attenuated signal at larger distances from the point of closest approach as illustrated in Figure[T] New 
hits on a given detector line are recursively sought within a narrow time window around the times of already selected 
hits. The expected hit times in adjacent or next-to-adjacent storeys are calculated, assuming the hits arrive linearly 
in the z-t plane as indicated by the asymptotes of the hyperbola in Figure [1] The projected arrival time f,+y at floor 
/ + j, (j - 1 , 2) is related to any pair out of the three hit times f,-, f/ipi , f,q:2 at floors /, / + 1 , / + 2 in the following way 

ti±j = ti + j(ti-ti^i), (2) 
ti±j = ti + j{ti-ti^2)/2, (3) 
ti±j = ffti + (; + 1) (f/Ti - fft2) ■ (4) 

In the following, time windows [f^"y'\fJ±J] based on the above equations are defined, which are used for the hit 
selection. If only two out of the three hits in floors / + 1, / + 2 exist, only one of the extrapolations above can be 
performed. A new hit in floor / + j is accepted if it occurred within a time window of 

= t^^J-P^ (5) 



f'f ' = max(fi -I- jAz- + f„ + jts), (6) 



n 
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Time (ns) Time (ns) 

Figure 1: The pattern of Cherenkov light from a muon track; (left: upward-going; right: downward-going) on a single detector line. Vertical 
position (z-coordinate) along the detector line and photon aixival time are related by a hyperbola (see Equations 117118) . The point z = 0,t = 
defines the point of closest approach between track and detector Hne; here, the track passes at 10 m. The colours correspond to different track 
angles: black 90° (vertical), red 60°, green 30°, blue 0° (horizontal). The thickness of the lines indicates the brightness of the expected signal. 



with f j and Az from Equation [T] This time window is asymmetric according to the generic shape of the Cherenkov 
cone hyperbola in the z-t plane (see Figure[T]and Equations I 171 1 SI discussed later). If all three hits exist, the minimum 
and maximum f,±j are determined and a time window of 

Cf = (7) 
= t'!:f + jts, (8) 

is used for floor ; + j. If several hits in floor / ± j meet the criteria, only the earliest one is selected. This procedure is 
applied recursively. If new hits cannot be found in either floor / + 1 or / + 2, the procedure stops, avoiding gaps of more 
than one floor The occurrence of several well-separated "hot spots" on the same detector line serves as independent 
seeds to restart the procedure. The described hit selection is performed once. It yields at most one hit per storey. The 
result of the hit selection on a bright neutrino candidate event from 2008 can be seen in Figure |3] 



6. Fitting procedure 

Before starting any fit, the list of selected hits is examined. Only events with more than 5 hits are accepted. If all 
selected hits are on a single detector line, a procedure for a single-line fit is started, otherwise a multi-line fit procedure 
is performed. Both procedures are discussed in the following. 

6.1. Fitting particle tracks 

A particle track is considered to be a straight line in space. The particle is assumed to move with the speed of light 
in vacuum. All space-time points, p(t), that are part of the track can be parameterized as 

pit) ^ q + c(t - to)u. (9) 

The particle passes through the point q at time to (i.e. q = p(to)) and moves in the direction u. The vector q can be 
shifted along the track by redefining to. Therefore the track is defined by a total of 5 parameters: three values to fix q 
for a given time and two angles to define m: 

M = {cos 0COS 0, cos 6?sin0, sinfl) , (10) 
5 



Cherenkov cone 




Track 



Detector line 

Figure 2: Illustration of the variables used to describe a track and its corresponding Cherenkov cone with respect to a vertical detector line. 

with 9 the elevation angle and (p the azimuth angle. The detector lines are approximated as vertical lines along the 
z-axis (see Section 2) at fixed horizontal positions (Lv, L,,). The z-component of the point of closest approach of a 
particle track to a detector line is given by 



Zc = 



q, — u,{q ■ it) + u,(LxUjc + LyW,.) 

1 -M? ' 



(11) 



through which the particle passes at time 



at a distance 



fc = fo + -(LxUx + LyUy + ZcUz - q - it), 
c 



dc = ^|{pAtc) - L,)^ + (Py(tc) - A.)2. 



(12) 



(13) 



The above-defined variables are illustrated in Figure |2] If the track is exactly vertical, and therefore parallel to the 
detector line, tc = fo and Zc - qz are chosen. The track can be conveniently reparametrized in terms of Zc, tc, dc for each 
detector line and the two angles which define tt. For a single-line fit the detector line can be placed at the coordinate 
origin (Li, Ly) = (0, 0) and Equations ll Ill2ll3] simplifv to 



tc - to + 



q- - Uziq- it) 

1 -M? ' 
1:Mz ~ (4 ' ^ 



c(l - M?) 



dc = ■y[p]<Jc)Tp^). 



(14) 

(15) 
(16) 
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A single detector line is invariant against a rotation of the coordinate system around the z-axis. This means the track is 
not fully defined and only 4 parameters Zc, tc, dc and u. can be determined. No dependence from azimuth is observed. 

To perform the track fit, it is necessary to know (a) the arrival time ty of a Cherenkov photon at the detector line 
positions (Lj, L^-, z), (b) its corresponding travel path dy and (c) its inclination with respect to the detector line, cos 9y. 
All three values can be derived from the parameters defined above and the refractive index n which is related to the 
Cherenkov angle Qc by 1 /n = cos Oc'. 

tyiz) = (tc + \ ((z - + """"^rte)) ' (17) 



dy(z) = , " Jd'- + {z-Zc)Hl-ul), (18) 

cos0^(z) = (l-M2)iL££ + ^. (19) 
" dy{z) n 

These equations hold exactly for Cherenkov photons of a given wavelength. Dispersion and group velocity effects as 
well as delays due to light scattering are ignored. An appropriate effective refractive index is used, depending on the 
medium in which the detector is installed. 



6.2. Fitting point-like light sources: the bright- point fit 

In contrast to a particle track, a 'bright-point' is defined as a point-like light source which emits a single light flash 
at a given moment. The light emission is assumed to be isotropic. The model of a bright-point does not only apply to 
artificial light sources such as optical beacons installed for detector calibrations fil5il but also in some approximation 
to light from hadronic and electromagnetic showers, for which the actual extension of the shower is typically smaller 
than the detector line distances. A bright-point is defined by four parameters: its position q and its time fo. In analogy 
with the definitions of the point of closest approach for particle tracks, it is straightforward to see that for a bright-point 
Zc - 9:, tc - to and 

dc = yjiq, - U)^ + {qy - Ly)\ (20) 

which simplifies to 

dc = ^ql + qj, (21) 

for a single detector line at (L^, Ly) = (0, 0). In this case, only the three parameters Zc, tc, dc can be determined, which 
again means that the number of parameters is reduced by one. The photon arrival time ty, its travel length dy and its 
angle with respect to a given arrival point z along the detector line can thus be determined in analogy to the case of a 
particle track: 

dyiz) = ^dl + {z-qd\ (22) 

tyiz) = fo + -dy, (23) 
c 

cos0.(z) = (24) 

dy 

As in the treatment of particle tracks, it is assumed that all photons have the same wavelength, which corresponds to 
a fixed effective refractive index. 



6.3. Quality function 

The quality function is based on the time differences between the hit times f, and the expected arrival time ty 
of photons from the track or bright-point, as in a standard x~ fit- The quality function is extended with a term that 
accounts for measured hit charges a,- and the calculated photon travel distances dy. The full quality function is 



Ni,i, 



1=1 



(ty - tif A(ai)D(dy) 
— 1 — 

en {a)dQ 



(25) 



where (a) is the average hit charge calculated from all hits which have been selected for the fit. The normalisation 
term do and the functions A(fl,) and D{dy) are discussed below. The second term of the right hand side of equationlZSl 
exploits the fact that an accumulation of storeys with high charges (hot spots) is expected on each detector line at its 
point of closest approach to the track. If such hot spots on several detector lines are arranged in a way that their z, t 
coordinates indicate an upward-going pattern, the event has indeed a high probability to originate from an upward- 
going neutrino. This allows the isolation of a high purity neutrino sample. The concept is illustrated in Figure |3]where 
the arrangement of the "hot spots" on the various detector lines allows to classify this event as upward-going even 
without an overlaid fitted track. The quality function depends on a number of parameters which have been tuned for 



sea water |112|,[lJ and the ANTARES detector ffl. 

The timing uncertainty, cr, is set to 10 ns for a, > 2.5 photoelectrons and to 20 ns otherwise. These values are 
significantly larger than the intrinsic transit time spread of the ANTARES phototubes (1.3 ns). They take into account 
an additional smearing due to the applied geometrical approximations from Section 3 and due to the occasional 
presence of late hits from small angle scattering or from electromagnetic processes which passed the hit selection. 

The charge term of the quality function is not arranged as a difference between expected value and measured 
charge in order to avoid incorrect penalties from large expected charges. Instead, the chosen form gives a heuristic 
penalty to the combination of high charge and large distance. The product is divided by the average charge of all the 
selected hits in an event, (a), to compensate for the fact that more energetic tracks or showers will produce more light 
at the same distance. The normalisation do balances the weight between the two terms. For ANTARES a value of 
c/o = 50 m has been chosen, motivated by the fact that at this distance the typical signal in a photodetection unit which 
points straight into the Cherenkov light front is of the order of one photoelectron. 

The charges have to be corrected for the angular acceptance of the storey. Assuming a downward-looking hemi- 
spherical photodetection surface with a homogeneous efficiency for photodetection (a good approximation for an 
ANTARES storey), the following simple correction function for the hit charges a', is obtained: 

' (26) 



cos 0y + I 



The average charge {a} is calculated from these corrected hit charges: 

<a) = y a'. (27) 

Furthermore, charges are protected against extreme values by the following function: 

aoa'. 

Aiad = , ' ■ (28) 



The function A(a,) introduces an artificial saturation such that for a| <s: oq, the charge is relatively unaffected, i.e. 
A(a,) X fl', whereas for a'. » oq, the charge saturates at A(fl,) x qq, the chosen saturation value. 
The photon travel distance is similarly protected using the function 



D(dy) = ^dj + d^, (29) 



which introduces a minimal distance d\. For large distances, dy » c/i, D(dy) ~ dy, whereas for very small distances, 
dy ^ di, D{dy) ~ di. This avoids an excessive pull of the fitted trajectory towards the detector line. 

Equations|28]and|29]can be motivated by the following consideration. The intensity of Cherenkov light decreases 
linearly with distance (ignoring absorption, dispersion and similar effects). Thus one can write a'^dy - aodi. For 
ANTARES, the constant a^di should be around 50 mxphotoelectron, which corresponds to the observation that about 
50 photoelectrons can be measured for a minimum ionizing particle at 1 m distance, or equivalently 1 photoelectron is 
seen at a distance of 50 m. Using this identity, it follows that A(fl,)D(c/y) = a'.dy, i.e. the two functions have no effect 
on the product of charge times distance for Cherenkov light; however, for background light they produce the desired 
penalty. 
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6.4. Minimization procedure 

The MIGRAD function of the MINUIT package fl^ is used to determine the minimum of the quaHty function Q. 
The value of Q/NDF (with NDF the number of degrees of freedom of the fit) after minimisation is called^? quality 
Q and is used for further selections. For the single-line bright-point fit, the parameters determined are the point of 
closest approach {i.e. Zc, tc, dc) with an additional parameter u^ in the case of a track fit. Equations [T7lfT9l and I22H24] 
are used to obtain dy, ty and cos 9y for particle tracks and bright-point fits, respectively. 

The multi-line fits use the same set of equations. But in this case one more parameter is needed to evaluate them. 
For the multi-line track fit q, m- and (p are used, for the multi-line bright-point fit q and tc- Table [1] summarizes the 
number of fit parameters for the different cases. 





single-line 


multi-line 


particle track 
bright-point 


4 
3 


5 
4 



Table 1 : Number of parameters for the different fits. Since multi-line events do not have rotational symmetry, an additional parameter is required. 



A critical issue is the choice of the starting values and possible boundary conditions. As explained above, the set 
of fit parameters used for single-line bright-point fits is Zc, tc, dc- Considering that most of the hits are grouped around 
the point of closest approach, Zc ~ (z/) at t^ ~ {ti), the mean values of the selected hit positions and times are used 
as starting values. The initial choice of dc is not critical for the fit stability. It is set to dc - 10 m. This completely 
determines the starting values for the single-line bright-point fit. For the single-line track-fit, a starting value for u, 
has to be chosen as well. To avoid finding a wrong u~, due to the inherent ambiguities when fitting Cherenkov light 
cones, three different starting values of M; are tried; -0.9, and 0.9. The fit with the smallest fit quality Q is taken. 
The allowed parameter range for u^ is limited to physical values {i.e. -1 < M; < 1). Fit results for which the fit stops 
at a boundary are excluded. 

The multi-line bright-point fit starts from the center of gravity of the selected hits in space and time. The multi-Une 
track fit requires, however, a prefit. For this purpose, an improved version iH of the "DUMAND prefit" 1 17| is used, 
which fits a straight line to the hit positions while allowing for a variable effective particle speed. 



7. Test of the algorithm with Monte Carlo simulations and ANTARES data 

The entire algorithm as described in Sections 3-6 was tested with ANTARES data and simulations of atmospheric 
muons and neutrinos. 



7.1. Data sample 

For the following analysis data are used which have been taken from December 2007 until end of 2008. During 
this time a minimum of 9 detector lines was active. From May 2008 on the ANTARES detector was complete with 
its 12 lines. Only runs are used, which have at least 80% of active channels when averaged over the run. A channel 
is called active if it has an instantaneous counting rate of at least 40 kHz, well below the typical baseline noise rate of 
60 kHz. Further, runs are excluded for which more than 40% of the active channels have been in burst mode, which 
means they measure an instantaneous counting rate more than 20% above the usual baseline rate. This selection 
provides a data sample which corresponds to an active livetime of 173 days. 



7.2. Simulations 



Downward-going atmospheric muons were simulated with Corsika 11 1811 . The primary flux was composed of 
several nuclei according to ( 19| and the QGSJET hadronic model ll20ll was used in the shower development. An 
average Mediterranean atmosphere is simulated, i.e. one which corresponds to a situation typically found in spring 
or autumn. Seasonal variations of the muon flux due to atmospheric temperature and density variations are found to 
be smaller than 3% in this region 112 ill , and they are ignored in this analysis. Upward-going neutrinos were simulated 
according to the parameterization of the atmospheric flux from l22ll in the energy range from 10 GeV to 10 PeV. 
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Only charged current interactions of and were considered. Disappearance of muon neutrinos due to neutrino 
oscillations was included in a simplified two-flavor model assuming maximal mixing and Am- = 2.4 ■ lO^-'eV-. 

The Cherenkov light, produced in the vicinity of the detector, was propagated taking into account light absorption 
and scattering in sea water lll3ll . The angular acceptance, quantum efliciency and other characteristics of the PMTs 
were taken from and the overall geometry corresponded to the layout of the ANTARES detector f^. The optical 
noise has been variable during the considered data taking period. It is simulated from counting rates observed in 
real data. At the same time the definition of active and inactive channels has been applied from real data runs of the 
selected period. The generated statistics corresponds to an equivalent observation time of 10 years for atmospheric 
neutrinos and, depending on primary cosmic ray energy, from 2 weeks up to 1 year for atmospheric muons. 

A dedicated simulation has been done to study the impact of coincident atmospheric muon events. One month 
of equivalent livetime has been simulated using the MUPAGE code ||23ll . An increase of 0.1% is found at trigger 
level. The fraction of downward reconstructed tracks increases by a similar fraction whereas for misreconstructed 
upward-going tracks an increase of 0.2% is found without using a cut in the fit quality. 



7.3. Results 

The accuracy in the reconstruction of the elevation angle (the angle with respect to the horizontal plane), which 
needs to be good for up-down separation of the reconstructed tracks, is illustrated in Figure |4] The reconstructed 
elevation angle is plotted versus the true elevation angle for down-going atmospheric muons and upward-going at- 
mospheric neutrinos, that are reconstructed as multi-line events. No further event selection was applied. More than 
95% of the fits converge for all events with selected hits on at least two lines. Most events are located within a narrow 
band around the diagonal and 80% of the events have their elevation angle reconstructed to better than 5°. No other 
structures are visible on the plots. It should be noted that this excellent agreement is obtained, even though more than 
half of the triggered atmospheric muon events are muon bundles, whereas here just a single track reconstruction is 
performed. 

Figure |5] confronts data and the two Monte Carlo samples for the full 173 days livetime. The left plot of Figure |5] 
shows the fit quality for the multi-line track fit of all the tracks reconstructed as upward-going. A cut on the track 
fit quality can effectively suppress background from misreconstructed atmospheric muons while maintaining a good 
efficiency for the atmospheric neutrino sample. A cut in the track fit quality of Q < 1.4 ensures a 90% purity while 
keeping 48% of the total sample of upward reconstructed multi-line neutrino tracks. The cut selects 665 neutrino 
candidates in data (about 4 events per day), to be compared to 609 events from the atmopheric neutrino Monte Carlo 
sample and 40 misreconstructed downward-going atmospheric muon events. No event from the coincident muon 
sample survives the cut which gives a 90% upper limit of 13 such events in the final sample. Data and simulation 
agree within the shown statistical errors of the data for Q < 1 .4, whereas for Q > 1.4 there is a systematic 30% excess 
of data. This excess is well within the estimated systematic error of the atmospheric muon simulation (see below). 
The resulting elevation angle distribution of the selected events with Q < 1 .4 is shown on the right plot of Figure [5| 
An excellent agreement between data and the atmospheric neutrino Monte Carlo sample based on the Bartol flux ||22|| 
is observed for the upward-going events, whereas the downward-going part features again a 30% excess of data with 
respect to the atmospheric muon simulation. The grey band gives the systematic error of the simulation. An extensive 
discussion of systematic effects is given in Is^l- Main contributions are from the uncertainty in the effective area of the 
PMTs (10%), the uncertainty in the water absorption length (10%) and the uncertainty in the angular acceptance of 
the PMTs. The latter is better known close to the PMT axis (15% error) than in the backward hemisphere (up to 35% 
error) which leads to a larger total systematic error of 45% for downward-going atmospheric muons than for upward- 
going atmospheric neutrinos (32%). Within the quoted errors, data are compatible with the chosen flux models both 
for atmospheric neutrinos and for atmospheric muons. 

Figure |6] compares data and simulation of some related quantities as the number of storeys used in the fit and the 
total charge of all hits used in the fit. Both plots show the comparison for the selected sample Q < 1 .4, which is 
dominated by atmospheric neutrinos. The shape of both distributions in data is well reproduced in the simulations. 

Figure|7]shows the effective area for all reconstructed multi-line events as function of the true neutrino energy. The 
comparison of the solid and the dotted line indicates that more than 40% of the events pass the quality cut of Q < 1 .4 
for energies below 10 TeV. This fraction reduces to 20% at 100 TeV and decreases further for higher energies. The 
low efficiency at highest energy can be explained by the use of a quality cut, which has been derived from Figure|5]for 
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Figure 3: Event display of a bright neutrino event from the 2008 data taking period. The arrangement of the 12 histograms corresponds approxi- 
mately to the octagonal layout of the 12 detector lines on the sea floor For each line the vertical position above the sea floor (m) is given as function 
of the hit time (ns) (see also Figure[T). r = corresponds to the time of the first hit which participated in the trigger. A window of -1000/ + 2000 ns 
is shown with respect to this time. ; = indicates the sea floor Active detector elements are between z = 100 and z = 450 m. Circles denote hits 
which participated in the trigger, crosses are other hits, boxes indicate the hit selection of Section[5] Colors refer to the hit charge as shown in the 
legend on the top light part of the figure. The line width and style of the fitted track illustrates the minimal distance between track and detector line; 
thick and solid lines stand for closer distances than thin or dotted fines. 
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Figure 4: CoiTelatioii between reconstructed and simulated elevation angle; Left: down-going atmospheiic muon events; Right: upward-going 
atmospheric neutrino events. 




Figure 5: Left: track fit quality for all upward reconstructed multi-line tracks for 2008 data (points with error bars) compared to a Monte Carlo 
sample of downward-going atmospheric muons (dashed histogram) and upward-going atmospheric neutrinos (sohd histogram); right: Elevation 
angle distiibution for events with Q < 1.4. The grey band indicates the error band of the combined prediction from neutrino and muon simulations. 




Figure 6: Left: Number of storeys and Right: total charge used in the track fit for all upward reconstructed multi-hne tracks for events with 
Q < 1 .4. Data (points with eiTor bars) are compared to a Monte Carlo sample of downward-going atmospheric muons (dashed histogram) and 
upward-going atmospheric neutrinos (solid histogram). 
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Figure 7: Effective area for all multi-line events (solid), selected events with Q < 1.4 (dotted) and selected events with Q < 1.3 + (0.04 ■ NDF)- 
(dashed) as function of the true neutrino energy for upward- going neutrinos. 

atmospheric neutrinos. By modifying the selection criteria one can easily recover the high energy part. The condition 
Q < \3 + (0.04 ■ NDF)^ is used for the dashed line. It has the same performance as 2 < 1.4 for Ey < 10 TeV 
where atmospheric neutrinos dominate, but a much improved performance at > 10 TeV. NDF - Nswrey - 5 is 
here used as a simple estimator for the brightness and therefore energy of an event. This modified cut selects 624 
events in data to be compared to 588 atmopheric neutrinos and 40 atmospheric muons keeping the quoted 90% purity 
of the neutrino sample. A full likelihood fit |Q] has a similar performance at high energies when requiring a 90% 
purity for atmospheric neutrinos. For energies below 1 TeV the method presented here is found more efficient than 
the likelihood fit. 




Figure 8: The distributions, normalized to the same content, of the angular error in degrees for selected atmospheric neutrinos (full) and atmospheric 
muons (hatched); left: elevation angle eiTor for multi-line events; right; space angle eiTor for non-coplanar hit selection. 

Figure |8] shows the angular resolution of the data sample which is selected by a cut of 1 .4 on the fit quality. 
Both plots of Figure[8]are integrated over the full angular range and over the atmospheric neutrino energy spectrum. 
The elevation angle is reconstructed with a precision of 0.4° for upward-going neutrino-induced muons and 0.5° for 
down-going muons (median of the angular error compared to the true muon direction). The space angle can only be 
unambiguously defined if the selected hits do not lie in a plane. Mirror solutions in azimuth are otherwise unavoidable. 
Because of the simplifications of the detector geometry used during reconstruction, a non-coplanar hit selection can 
only be achieved when hits from at least three detector lines are used. This condition is fulfilled for 25% of the 
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selected multi-line atmospheric neutrino events. For this sub-sample a median space angle error of 1.0° is found for 
atmospheric neutrinos and 1.4° for atmospheric muons as shown on the right plot of Figure |8] 



7.4. Improving the angular resolution 

The angular resolution which has been quoted in the previous section can be improved in several ways, without 
changing the event selection and without compromizing the speed of the algorithm. 

Resolving azimuthal degeneracy 

For a particle trajectory reconstructed using hits on only two straight detector lines, there always exists an alterna- 
tive trajectory having an identical g-value, but a different direction. The two trajectories will have the same elevations, 
but will differ in their azimuthal orientation. The degenerate trajectory is easily determined given the positions of the 
two detector lines and the trajectory resulting from the fit. First, the point at which the particle trajectory crosses the 
plane connecting the two detector lines is determined. The trajectory is then rotated about this point by an angle lif/ 
to form the degenerate track candidate, where if/ is the angle between the fit result and the plane connecting the two 
detector lines. 

Armed with both possible track solutions, the two candidates can be discriminated with the following procedure. 
The time residuals of all hits are calculated from the differences of the actual hit times and the theoretical hit times 
according to Equation [17] The sum of the charge of each hit having a time residual smaller than 30 ns is then 
calculated. This can include hits from any line and not just the two used in the original fit, helping to break the 
degeneracy. The hit selection is done separately for both track candidates and the candidate with the largest total 
charge is chosen. In simulations of both downward-going muons and upward-going neutrinos this procedure chooses 
the correct track candidate 70% of the time and improves the two-line angular resolution by 42% for upward-going 
atmospheric neutrinos and by 15% for downward-going atmospheric muons. 

Additional fit step 

Taking the result of a multi-line track fit or its mirror solution as it was described in the previous sections as a 
prefit, a new hit selection and a subsequent fit is performed. Every hit which has a time residual smaller than 20 ns 
with respect to the prefit track is chosen. These hits are then used to minimize the function 




{ty - tjf 
2cr2 



(30) 



which is a robust estimator that combines the properties of both least-squares and absolute- value minimizers 02411 . 
This estimator is quadratic for small values of the time residual {i.e. when \ty — t\jcr <K 1), just as with a standard 
estimator, but is linear for large time residuals. It is therefore not affected as dramatically by background hits with 
large residuals surviving the hit selection. Note that only the time residual of hits are minimized at this stage; no 
information concerning the number of photoelectrons detected is used. For the current studies cr = 1 ns has been 
taken. The precise value of cr has little impact on the angular resolution. 



Use of the detailed detector geometry 

While the algorithm has so far assumed a simplified detector geometry, it is possible to extend it in order to 
exploit the actual location and orientation of the PMTs. This detailed geometry information can be employed in the 
calculation of the predicted photon arrival time {ty in Equation |30]|. The use of the resulting time residual for both the 
hit selection and track fit minimization, as described in the previous section, has been studied. Simulations for this 
study have been performed in the absence of sea currents, i.e. again for straight detector lines. The systematic angular 
shift due to the coherent inclination of all detector lines in a horizontal sea current is estimated to be smaller than 0.2° 
for a typical sea current of 5 cm/s as measured at the ANTARES site. However, such a systematic shift is not included 
in the numbers and figures quoted below. 
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Results 

The effect of these three additional steps can be seen in Table |2] which compares the median angular error of the 
multi-line algorithm described in Section ISTI to the subsequent fit using Equation [30l with and without the use of a 
detailed detector geometry. The condition on the number of detector lines A^, from which hits have been used in the 
fit, is specified in the first column of Table |2] Results are given for upward-going atmospheric neutrinos as well as 
for neutrinos from an E:^~ flux, which is typical for astrophysical sources. Only events which fulfill the condition 
e < 1.3 + (0.04 ■ NDF)~ are included. 

For tracks reconstructed using hits on only two lines in the initial fit, the improvement can be attributed mainly to 
the partial resolution of the azimuthal symmetry. For the case where more than two lines had been used initially, the 
main improvement comes from the additional fit step. An additional 0.1° in resolution may be gained by using the 
detailed detector geometry within each storey. The fuU multi-line sample (A^ > 1) profits from a combination of both 
effects. 





a for atmospheric v 




a for 


oc E-'^ 


lines 


Eq.|25j 


Eq.l3()| 


Eq.l30l(geom) 


Eq.|25] 


Eq.[30] 


Eq.l30l(geom) 


N=2 


4.2° 


1.7° 


1.5° 


4.4° 


1.3° 


1.1° 


N>1 


2.9° 


1.3° 


1.1° 


1.3° 


0.60° 


0.52° 


N>2 


1.00° 


0.64° 


0.54° 


0.94° 


0.47° 


0.41° 



Table 2: The median angular error (a) in the multi-line fit algorithm as of Section ISTI and the subsequent fit with and without using geometry 
information, for atmospheric neutrinos and neutrinos from an E^- flux. 

The improvement in the angular error for neutrino-induced muons can also be seen in Figure|9](left), which shows 
the energy dependence of the median angular deviation from the true muon direction for all reconstructed tracks 
using hits on more than 2 lines having a g-value better than 1.3-1- (0.04 ■ NDF)^. The selection is applied only to 
the original track fit, so that the observed improvement of the angular resolution is the result a refined hit selection, 
combined with a new fit. The three curves correspond to the multi-line track fit from Section |6T| (solid) and the fits 
using Equation[30]with an approximated geometry (dotted) and with a detailed geometry of each storey (dashed). The 
muon track angular error is almost energy independent at neutrino energies above 10 TeV. It converges to the values 
quoted in table|2]for an E^^ flux. 

This can be compared to the performance of a full likelihood fit |@] which uses the full geometry, a detailed 
modelling of the time residual, {e.g. an energy dependent tail of late hits and a contribution of noise hits) and which 
applies a scan of starting positions of the fit to avoid local minima. Such a procedure can reach an angular resolution 
of 0.3° at high energies for the entire multi-line sample, while using about 10 times more computing time than the 
method presented here. 

Figure |9] (right) shows the point spread function for neutrino events after the fit from Equation|30] Contrary to the 
figures shown earlier the angular error is given here with respect to the neutrino direction - a quantity which is relevant 
to point back to a potential neutrino source. The optical followup program ^ uses a sub-sample of the events, shown 
on the dashed line in Figure |9] (right). As the field of view of the used robotic telescope is of the order of 2°, it can be 
seen that more than 70% of the candidate sources would lay inside the field of view. 

The improvement is further demonstrated by the time residual distributions, as shown in Figure [10] The time 
residuals are shown for all hits in each event for atmospheric neutrino events. The residuals of the subsequent fit using 
the detailed geometry show a more pronounced peak at zero, whose width is roughly 40% smaller compared to the 
multi-line fit from Section ISTI 

These extensions to the algorithm do not significantly affect the time needed for track reconstruction. Including 
the subsequent fit, the entire reconstruction algorithm takes an average of 5 ms per muon event on an Intel Xeon 
2.66 GHz CPU for events recorded by the 12-line ANTARES detector. 

8. Conclusions 

A generic and fast track fit algorithm for muon tracks in neutrino telescopes has been presented which can reliably 
distinguish upward-going neutrinos from the overwhelming background of downward-going muons. The method is 
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Figure 9: Left: Dependence of the angular resolution for selected tracks (Q < 1.3 + (0.04 ■ NDF)^,N > 2) on E,.. The median angulai' en'or 
of the multi-line fit (solid line) is compared against the fit from Equation [30] without (dotted) and with (dashed) use of a detailed description of 
the geometry of a storey; right: Cumulative point spread function for selected upward-going neutrino events with an E^^ spectrum and Q < 
1.3 + (0.04 ■ NDF)",N > 2 for multi-line fits (solid) and fits from Equation l30l with an approximate geometry (dotted). 

used in ANTARES for different applications such as an online event display, a neutrino monitor and an alert sending 
program to trigger optical follow-ups of selected neutrino events. It has also been used for an analysis of atmospheric 
muons Js]] ■ Other potential applications include the study of magnetic monopoles and nuclearites and the measurement 
of atmospheric neutrinos. 
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